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A  PERFECTLY  MATCHED  LAYER 
FOR  PERIDYNAMICS  IN  TWO  DIMENSIONS 

Raymond  A.  Wildman  and  George  A.  Gazonas 

A  perfectly  matched  layer  (PML)  absorbing  boundary  is  formulated  for  and  numerically  applied  to  peri- 
dynamics  in  two  dimensions.  Peridynamics  is  a  nonlocal  method,  derived  to  be  insensitive  to  discontinu¬ 
ities,  more  easily  simulating  fracture.  A  PML  is  an  absorbing  boundary  layer,  which  decays  impinging 
waves  exponentially  without  introducing  reflections  at  the  boundary  between  the  computational  region 
and  the  absorbing  layer.  Here,  we  use  state-based  peridynamics  as  PMLs  are  essentially  anisotropic 
absorbing  materials,  therefore  requiring  arbitrary  material  parameters.  State-based  peridynamics  is  also 
more  convenient  for  auxiliary  field  formulations,  facilitating  the  implementation  of  the  PML.  Results 
show  the  efficacy  of  the  approach. 


1.  Introduction 

Originally  introduced  in  [Silling  2000],  peridynamics  is  a  nonlocal  formulation  of  elastodynamics,  which 
can  more  easily  incorporate  discontinuities  such  as  cracks  and  damage.  Derivatives  of  field  variables 
in  the  classical  continuum  model  are  replaced  by  integrals  over  a  small  neighborhood  of  microelastic 
kernels  that  replace  standard  constitutive  relations.  In  its  discretized  form,  an  elastic  solid  is  treated  as 
a  collection  of  particles  or  nodes,  each  connected  to  its  neighbors  by  breakable  bonds.  Bond  breakage 
can  be  defined  to  occur  when  a  bond  is  stretched  past  some  predetermined  limit.  After  a  bond  is  broken, 
any  supported  force  transfers  to  the  remaining  bonds,  increasing  their  supported  load,  and  encouraging 
more  breakage.  Eventually,  this  process  autonomously  leads  to  cracking  and  failure.  The  end  result  is 
a  method  capable  of  predicting  crack  growth  in  brittle  elastic  materials  [Gerstle  et  al.  2005;  Silling  and 
Askari  2005;  Emmrich  and  Weckner  2006;  Demmie  and  Silling  2007;  Kilic  et  al.  2009;  Ha  and  Bobaru 
2010], 

Over  the  last  decade,  peridynamics  has  been  extended  past  its  original  formulation.  First,  the  nu¬ 
merical  method  originally  outlined  in  [Silling  and  Askari  2005]  has  been  extended  to  include  adaptive 
refinement  [  Bobaru  et  al.  2009],  replaced  with  different  quadrature  rules  [Emmrich  and  Weckner  2007], 
and  implemented  in  a  parallel,  molecular  dynamics  code  [Parks  et  al.  2008].  In  addition,  it  has  been 
extended  to  different  material  types  including  viscoplastic  [Foster  et  al.  2010],  micropolar  [Gerstle  et  al. 
2011],  and  nanofiber  networks  [Bobaru  2007].  It  has  also  been  applied  to  different  fields  such  as  heat  con¬ 
duction  [Bobaru  and  Duangpanya  2010]  and  electromigration  [Gerstle  et  al.  2008].  Aside  from  practical 
applications,  the  mathematics  behind  the  approach  have  been  studied:  Weckner  et  al.  [2009]  derived  a 
Green’s  function  for  the  peridynamic  equation  and  Weckner  and  Abeyaratne  [2005]  discussed  dispersion 
relations  for  various  kernels.  Most  importantly  for  this  work,  state -based  peridynamics  was  introduced, 

Keywords:  peridynamics,  perfectly  matched  layer,  absorbing  boundary. 
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allowing  for  more  flexible  constitutive  relations  [Silling  et  al.  2007],  As  will  become  clear  later,  state- 
based  peridynamics  allows  for  an  auxiliary  field  formulation,  which  is  necessary  for  the  implementation 
of  a  perfectly  matched  layer  (PML). 

While  most  peridynamics  work  has  focused  on  simulating  problems  with  free  or  fixed  boundary 
conditions,  there  are  applications  in  which  the  simulation  of  an  infinite  medium  may  be  useful,  such 
as  wave  or  crack  propagation  in  a  half-space.  Absorbing  boundary  conditions  are  a  way  of  simulating 
an  infinite  medium  by  absorbing  any  impinging  waves  at  the  computational  boundaries  so  they  do  not 
reflect  back  into  the  simulation.  A  PML  is  such  an  absorbing  boundary,  and  was  originally  introduced 
for  electromagnetic  simulations  [Berenger  1994;  Chew  and  Weedon  1994].  PMLs  differ  from  traditional 
absorbing  boundary  conditions  in  that  they  are  an  absorbing  layer,  placed  between  the  computational 
region  of  interest  and  the  truncation  of  the  grid  or  mesh.  They  can  also  be  thought  of  as  an  anisotropic 
absorbing  material,  which  is  why  the  flexibility  of  a  state-based  peridynamics  is  necessary. 

PMLs  have  two  important  qualities:  First,  waves  in  a  PML  decay  exponentially,  and  second,  in  their 
analytic  form,  no  waves  reflect  at  the  interface  of  a  PML  and  the  computational  region.  These  properties 
make  them  ideal  for  simulating  wave  propagation  in  infinite,  unbounded  regions.  Since  their  introduction, 
PMLs  have  been  extended  to  many  different  types  of  media  [Uno  et  al.  1997;  Teixeira  and  Chew  1998; 
Dong  et  al.  2004],  different  numerical  methods  [Pissoort  and  Olyslager  2003;  Pissoort  et  al.  2005;  Alles 
and  van  Dongen  2009],  and  different  fields  [Chew  and  Liu  1996;  Liu  and  Tao  1997;  Festa  and  Nielsen 
2003]. 

This  paper  implements  a  peridynamic  formulation  of  elastodynamics  in  two  dimensions  and  terminates 
the  boundary  with  a  PML.  As  is  discussed,  the  use  of  a  PML  is  facilitated  with  an  auxiliary  field  formu¬ 
lation,  derived  from  state -based  peridynamics,  and  the  peridynamic  equation  is  broken  into  five  coupled 
equations.  A  PML  was  applied  to  one-dimensional  peridynamics  in  [Wildman  and  Gazonas  201 1],  which 
used  the  results  of  [Du  et  al.  2012]  to  formulate  an  auxiliary  field  equation.  This  approach  required  a 
matrix  representation  of  the  auxiliary  field,  which  may  be  memory  prohibitive  in  higher  dimensions. 

The  remainder  of  the  paper  is  organized  as  follows:  Section  2  discusses  the  formulation  of  peridynam¬ 
ics,  PMLs,  and  their  numerical  implementation;  Section  3  gives  some  results;  and  Section  4  summarizes 
the  report  and  details  future  work. 


2.  Formulation 


In  this  section,  a  PML  is  formulated  for  state -based  two-dimensional  peridynamics.  First,  in  Section  2A, 
a  linear  elastic,  state-based  peridynamics  formulation  will  be  reviewed.  Next,  Section  2B  reviews  the  for¬ 
mulation  of  a  PML.  Section  2C  then  applies  the  PML  to  state-based  peridynamics,  and  finally  Section  2D 
discusses  a  discretization  of  the  formulation  using  the  standard  node -based  peridynamics  method. 


2A.  Two-dimensional,  state-based  peridynamics.  The  continuum  equation  of  motion  in  an  elastic  solid 
can  be  stated  as 


•  a  +  b, 


(2-1) 


where  (in  two  dimensions)  p(x)  [kg/m2]  is  the  density,  u(x,  t )  [m]  is  the  displacement,  d(x.  t)  [N/m]  is 
the  stress  tensor,  and  b(x )  [N/m2]  is  a  body  force  [Malvern  1969].  (Throughout,  boldface  type  denotes 
a  vector  and  a  boldface  variable  with  an  overbar  denotes  a  tensor.)  Equation  (2-1)  is  a  local  formulation 
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because  the  divergence  of  the  stress  (and  gradient  of  the  displacement  implied  in  its  definition)  represents 
a  local  operation  on  a  variable.  In  other  words,  the  action  of  V  •  a  only  depends  on  a  at  a  single  spatial 
point.  In  problems  involving  discontinuities,  such  as  cracks,  the  divergence  at  such  discontinuities  is 
not  well  defined,  leading  to  numerical  implementation  problems.  Peridynamics  proposes  replacing  V  •  a 
with  a  nonlocal  operation  that  nonetheless  also  represents  a  force 


d2  f 

p — =  /  f(u  —  w,  x  —  x)dVx'  +  b, 
dt  Jur 


(2-2) 


where  fix'  —  x,  u'  —  u)  [N/m4]  represents  a  micromodulus  force  function  (or  kernel)  that  defines  a 
force  between  two  points  and  Htx  represents  a  horizon  or  maximum  distance  over  which  two  points  can 
influence  each  other  [Silling  2000],  The  micromodulus  function  becomes  the  constitutive  response  in 
the  formulation,  replacing  Hooke’s  law  in  the  continuum  case.  In  its  original  form,  the  micromodulus 
function  was  developed  as  a  simple  elastic  response  following 


/0?>$)  =  c 


h  +  H  \$  +  ri\-\y\ 
\%+v\  \v\ 


(2-3) 


where  c  is  some  constant.  Hi  ■ )  is  the  Heaviside  step  function,  and  S  is  the  radius  of  the  horizon  region, 
which  defines  [Silling  2000].  Note  that  in  contrast  to  [Silling  2000],  for  simplicity  there  is  no  history- 
dependent  failure  term  in  (2-3).  Equation  (2-3)  is  isotropic,  though  not  strictly  linear  in  terms  of  u. 
Linearizing  (2-3)  gives 


f(V,S)  =  C($)ii,  CV-cM  Cm  =  cHi8-\i;\), 

l?l 


(2-4) 


where  ®  denotes  an  outer  product  [Silling  2000].  The  function  C(|£|)  is  the  kernel  function,  typically 
taken  to  be  a  Heaviside  function.  Here,  we  will  also  use  a  Gaussian  kernel  function,  which  tends  to  give 
smoother  results  with  less  apparent  ringing  in  the  solution.  A  Gaussian  kernel  is  defined  as 


CGauss(lll)  =  Ce-(l*l/5)2. 


(2-5) 


For  the  Gaussian  kernel,  the  horizon  8  does  not  delineate  a  strict  bond  family  as  the  Heaviside  function, 
but  describes  the  decay  of  the  kernel.  To  determine  the  bond  family,  Mx,  a  small,  arbitrary  value  can  be 
chosen  as  a  cutoff  for  the  kernel.  Note  that  the  cutoff  has  a  large  impact  on  the  efficiency  of  the  method: 
too  small  a  cutoff  and  a  large  number  of  bonds  must  be  included  in  each  calculation.  Throughout,  we 
set  the  cutoff  to  10-6. 

A  PML  application  requires  an  auxiliary  field  formulation,  as  it  is  essentially  an  anisotropic  absorbing 
material,  if  a  nonphysical  one.  Consequently,  a  state-based  peridynamic  formulation  [Silling  et  al.  2007] 
is  necessary  to  implement  the  required  constitutive  relations  in  the  absorber.  State-based  peridynamics 
uses  a  family  of  bonds  to  determine  a  given  force  rather  than  a  single  bond  independently.  This  more 
general  approach  allows  for  inelastic  behavior  and  more  general  elastic  behavior,  and  is  governed  by 

92  f  _  _ 

p—u=  /  (J[jc,  t](x'  —  x)  —  T[x' ,  t](x  —  x'))dVx>  +b,  (2-6) 

_  9a  J'mx 

where  T[ x,  t\(x'  —x)  is  a  peridynamic  vector  state,  with  the  parameters  in  the  square  brackets  indicating 
variables  that  act  as  arguments  to  any  functions  referenced  in  the  vector  state  and  the  variables  in  the 
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angle  brackets  acting  as  arguments  to  the  vector  state  itself.  In  the  state-based  formulation  of  [Foster 
et  al.  2010],  the  deformation  gradient,  given  by 

F  =  / +  «V,  (2-7) 


can  be  approximated  as  a  vector  state  as 

F[x,t]  = 

where  §  =  jc'  —  jc.  A"  isa  shape  tensor  given  by 


f  C(\i;\)(Y[x,t]{i;)®i;)dVxr 
LJkx 


K\ 


K[x,t]=  f  C(|||)(|®$)</Vxs 

Jwx 

and  Y  is  a  deformation  vector  state  given  by 

y[x,/]<$>  =  *  +  $, 


(2-8) 


(2-9) 


(2-10) 


with  r]  =  u\x' ,  f]  —  u\x,  t]  [Foster  et  al.  2010]. 

The  deformation  gradient  can  now  be  substituted  into  Hooke’s  law  and  strain-displacement  relations, 
giving  a  stress  term  a  in  terms  of  u  in  plane  strain 


where 


32  _ 

p—u  =  V  a  =  V-  (c  :  e), 


e[x,  t\  —  \  (Vm  +  mV)  =  i  (F'fjt,  t]  +  F[x,  t]T  —  21), 


c  = 


(1  +  v)(l  —  2v) 


1  —  v  v  0 

A  T  2/x 

A 

0 

v  1  —  v  0 

= 

A 

A.  2fi 

0 

i 

o 

o 

1 

K> 

0 

0 

2/z_ 

(2-11) 

(2-12) 

(2-13) 


E  is  the  Young’s  modulus,  v  is  the  Poisson's  ratio,  and  A  and  ji  are  the  Lame  parameters.  Ultimately, 
the  peridynamic  vector  state  T  for  plane  strain  elasticity  is  given  by 


(2-14) 


2B.  Perfectly  matched  layer.  The  first  step  in  formulating  a  PML  is  to  construct  an  analytic  continuation 
to  the  complex  plane 

x  =  x  +  ig(x),  (2-15) 

where  g(x)  is  a  given  function  describing  the  deformation  [Johnson  2010].  This  mapping  has  the  effect 
of  transforming  traveling  waves  of  the  form  elkx,  where  k  =  co/c  is  the  wave  number,  into  evanescent 
waves  of  the  form  e'kxe~kg(xK  thus  attenuating  such  waves  in  the  PML  region. 

Substituting  x  into  the  above  equations  would  yield  a  viable  method,  though  one  that  requires  complex 
coordinates.  A  simpler  solution  is  to  change  variables  back  to  the  real  paid  x,  which  requires  a  relation 
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for  the  differential  quantities,  given  as 


dx  = 


\  ,  •  d 
l  +  ,Txs. 


dx: 


(2-16) 


part  ial  differential  quantities  are  used,  as  this  is  used  as  a  substitution  for  the  above  equations  involving 
functions  of  both  x  and  t  [Johnson  2010].  A  convenient  choice  for  g(x )  is 

4>(x) 


TxSM  = 


CO 


(2-17) 


because  the  1  /co  factor  creates  a  frequency-independent  attenuation  rate  in  dispersionless  materials  [John¬ 
son  2010].  (Peridynamic  formulations  are  not  dispersionless,  as  discussed  in  [Weckner  and  Abeyaratne 
2005],  though  this  standard  choice  is  used  for  simplicity.)  Finally,  the  substitution  that  must  be  made  for 
any  spatial  derivative  can  be  written 


_3_  _ 1 _ 3_ 

dx  1  +  i<p(x)/co  dx ' 


(2-18) 


Before  applying  a  PML  directly  to  the  peridynamic  equation,  (2-1)  will  be  treated  so  that  the  PML  ap¬ 
plication  to  peridynamics  will  be  clear.  It  is  convenient  to  convert  (2-1)  to  the  Laplace  domain,  assuming 
e~st  time  dependence,  giving 

ps2u  =  V  a,  (2-19) 


where  the  Laplace  transform  of  a  variable  is  indicated  by  ££{/}  =  f.  Next,  we  express  the  wave  equation 
as  two  coupled  first-order  partial  differential  equations,  the  first  in  u  and  the  second  in  scj/  —  a: 

psu  =  V-i/r,  si/r  =  c:€.  (2-20) 


Expanding  (2-20)  into  components  gives  five  coupled  equations: 


3  7  .3  7  ~  3  7  .3 

pm'  =  +  TA"  ps“>  =  TA'  +  3 y*” 

'  "7  =  +  -Ii ) ii  ,  -  ^  (7,,  v 3,  +  VI.'/.  =  ((  ''  ■  ) - 

Here,  we  will  make  the  substitution  given  in  (2-18)  for  all  spatial  derivatives,  written  as 

_3_  ^  s  3 
dx  s  +  4>(x)  9* 

in  the  Laplace  domain,  and  later  define  cf>x  and  <py  in  the  desired  absorbing  boundary  locations. 
(2-21)i  as  an  example,  we  get 


(2-21) 


(2-22) 

Using 


psux 


s  3 
s  +  cj)x  dx 


Vri  + 


5  3  7 

s+4>y  3  y  T 


p(s  +  cpx)(s  +  (py)iix  -  (s  +  0v)  —  +  (s  +  <t>x)  —  fT. 


dx 


dy 


(2-23) 


The  remaining  components  of  (2-21)  can  be  expanded  in  a  similar  way. 

Wherever  cj)  >  0,  u  and  a  will  exponentially  decay.  Before  discretization,  any  change  in  cp  will  not 
result  in  any  reflections,  so  the  region  of  interest  would  have  </>  =  0,  and  the  PML  region  could  have  a 
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discontinuity  in  applying  cp.  In  practice,  however,  numerical  reflections  can  result  from  discontinuous 
material  parameters  after  discretization,  so  it  is  better  to  use  a  smooth  transition  for  cj).  Here,  we  divide 
the  PML  region  into  two  parts,  one  in  which  cp  is  a  constant  value,  and  the  other  in  which  cj)  ramps  up 
to  that  constant  value  following  a  Gaussian  distribution.  An  example  is  shown  in  Figure  1,  with  the 
constant  region  set  to  0.1  m  and  the  Gaussian  region  0.2  m.  The  variance  of  the  distribution  is  set  so  that 
the  minimum  value  in  the  Gaussian  region  is  I0  6. 


2C.  Auxiliary  field  formulation  and  PML  application.  Peridynamics  is  not  typically  stated  in  terms  of 
Cartesian  components  as  in  (2-24),  but  we  can  expand  the  state -based  formulation  into  components  and 
match  terms  to  (2-24).  Following  this  approach  yields  a  viable  method  for  performing  PML  substitutions. 
First,  the  state -based  peridynamic  equations  (2-6)— (2- 14)  can  be  written  explicitly  as 


psux\x,  s]=  [  C(|£|)[(Vc[x,  *  I  k'xx  +  ^x\x,  s\k™)ljx  +  0M*',  s]kr™  +  \jrT[x',  s]k'™)$x\dVX' 

J'Kx 

+  £  cm[(irx[x,  s\kf;  +  irT[x,  s]k™)$y  +  (&[*',  s\kf™  +  fx[ x',  s]k'™)^dVx,, 

psuy\x,  s]=  f  C(\$\)[(ijrT[x,  s \ kff  +  xjry[x,  s \ kff)^x  +  (xpT[x',  s \ k'f™  +  xfry[x',  s]k'™)i;x]dVx> 

JMx 

+  £  Cm[WAx,  S\kf;  +  j,y[X,  S]k™)$y  +  OMx',  S\kf-  +  l}y\x',  s]k$”)$y]dVx', 
sijirx[x,s]  =  (A +  2  fl) 


f  C(|§|)(yx[x,  s]$xk™  +  Yx[x,  s]%yklyf)dVx’  -  1 

Umx 


+  A 


f  C(|||)(yy[X,  S]$xk™  +  Yy[X,  S]$yk™)dVX'  -  1 

U 


S\jfy  [X,  S]  =  X 


f  C(|f  |)(y,[x,  s]$xk™  +  YX[X,  s]tjyk™)dVxl  -  1 
lJxx 


+  (A  +  2/x) 


f  c(i5i)(yy[x,  +  ry[x,  .stfyk'";)dvx.  -  i 

_Jxx 

4r[x,  s]  =  n  f  Cm{Yx[x.  s]ixk™  +  Yx[x,  s]ijyk™)dVx, 

JMx 

+  M  £c(|*|)(y,[x,  s]ixk™  +  Yy[x,  sViyk™)dVx,,  (2-24) 


where 


K1 


'Kim 

Kxx 

rim 

Lryx 


k 

k 


inv" 

xy 

inv 

yy  - 


(2-25) 


Though  no  derivatives  appear  in  (2-24),  the  correspondence  of  each  term  to  those  in  (2-21)  is  apparent  — 
with  partial  derivatives  following  the  component  of  £  —  and  the  PML  substitutions  can  be  made.  For 
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Figure  1.  An  example  of  a  PML  across  one  dimension  using  a  Gaussian  ramp. 


example,  the  first  equation  in  (2-24)  can  be  rewritten  as 

P(s  +  4>x)(s  +  <Py)ux  =  (s  +  4>y)  I  C(\%\)(fax,s]k™ +  fa[x,s]k™)l;xdVx’ 

+  (s  +  <l>y )  f  CQ$\){fax',s]kr™  +  fa[x',s]k?y™)SxdVX' 

Jvex 

+  (s  +  <px)  f  Cm(fax,s]k™  +  fax,s]k™)$ydVX' 

+  (s  +</),)  f  CQ$\)(fax',  s]kf™  +  fax',  s]k'™)tydVx>, 
with  the  remaining  equations  following  similarly. 


(2-26) 


2D.  Discretization.  For  the  temporal  discretization,  forward  Euler  will  be  used  to  simplify  the  presen¬ 
tation  and  implementation.  Higher-order  temporal  discretizations  can  be  used,  though  they  lead  to  more 
terms  in  what  follows.  The  Laplace  domain  was  used  throughout  to  facilitate  the  temporal  discretization 
of  the  final  equations.  Because  differentiation  in  the  Laplace  domain  is  represented  by  s,  approximations 
to  s  can  be  directly  substituted  in  terms  of  the  z-transform.  This  technique  is  used  in  the  design  of  digital 
filters,  where  it  is  known  as  “filter  design  by  approximation  of  derivatives”  or  “the  bilinear  transforma¬ 
tion”  [Proakis  and  Manolakis  1996],  and  in  integral  equation  methods  where  it  is  known  as  “convolution 
quadrature”  [Lubich  1988a;  1988b],  or  “finite  difference  delay  modeling”  [Wang  et  al.  2008].  A  forward 
Euler  approximation  can  be  stated  as  the  substitution 


5  — > 


z.  —  1 
At  ’ 


(2-27) 


where  z  is  the  z-transform1  variable  representing  a  unit  advance  and  At  is  the  time  step  size.  Substitution 


*The  z-transform  is  defined  here  as  X (z)  =  Z{x[n]}  =  YlnL o  •*[«] z  "  for  causal  signals. 
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into  (2-26)  gives 

p(z-l+At(l>x)(z-l+At(l>y)Ux  =  (z-l+Atfy)  f  Cm(^x\x,z\k^+^Ax,z\k^)^xdVx- 

+  (z-\+At<py)  f  Cm{*xlx\zlk'™+Vdx\z]k'™)$xdVX' 

+  (z-l+At<f>x)  I  Cm{Vx[x,z]k™+'l>T[x,z]k™)i;ydVx, 

J  df€x 

+  (z-l+At<px)  f  Cm{*Ax\z]k'™+VT[x',z]k'™)SydVxr,  (2-28) 

where  capital  letters  indicate  the  z-transform  of  a  variable.  Expanding  the  quadratic  term  on  the  left-hand 
side  gives 

(z  - 1  +  A t<px)(z - 1  +  A t(f)y)  =  A t2<px<py  +  A t<pxz  -  A t(j>x  +  At<f>yZ  +  z2  —  z  —  At<f>y  -  z  +  1 

=  z2  +  (At<px  +  At(f>y  —  2)z  +  A t24>xcf>y  -  A t<px  -  At<py  +  1.  (2-29) 

Multiplying  by  z-2  and  rearranging  gives  an  update  equation  in  terms  of  z,  which  can  be  converted 
to  a  time-stepping  method  via  the  inverse  z-transform2  (assuming  vanishing  initial  conditions  and  an 
appropriate  region  of  convergence)  as 

ux[x,  1 ]  =  —(yx  +  Yy)ux[x,  l- 1]  -  YxYyUx{x,  l- 2] 

+  —  f  CQ$\){l,x[xA-l]k™  +  fT[xJ-l]k™)$xdVX' 

P  Jwx 

+  —  f  C(|||)(^[x' ,1-lJC  +  *Ax\l-Ak’™%dVX’ 

p  Jwx 

+  —Yy  f  Cm{*AxJ-2]k™  +  xl,Ax,l-2\k™)$xdVX' 

P  Jmx 

+  —Yy  f  Cm(rAx\l-2]k'™  +  xl,Ax\l-2]k'™)ixdVx' 

P  '  Jvex 

+  —  f  cm)^Axj-\w:;+^Ax,i-A^;%dvxl 

P  Jytx 

+  —  f  Cm(xl,Ax\  l-Wk’™  +  iff  Ax',  l-l]k'™)$ydVx, 

P  J  'M,  x 

+  —  Yx  f  Cm{*Ax,l-2]k%  +  1rAx,l-2\k™)l;ydVX' 

P  Jwx 

+  —  Yx  f  Cm{rAx\l-2]k'™  +  xlsAx\l-2]k'™)t;vdVx,,  (2-30) 

P  Jvex 

2The  only  necessary  property  is  the  delay:  x[n  —  fc]  -o-  z~k X(z) 
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where  1  is  the  time  step  number,  yx  =  A t<j)x  —  1,  and  yy  —  At<py  —  1.  A  stress  component  update  equation 
becomes,  for  example. 


%//x[x,  l ]  =  ~(yx  +  Yy)r//X[x,  l- 1]  -  YxYyfxix,  1-2] 

A t  (A  -f-  2/z) 


f  C(|?|)(y,[x,  l-mk™  +  YX[X,  l-mk™)dVx ,  -  1 
\.J-xx 

f  C(Hi)(yv[x,  i-2]$xk™  +  yx[x,  i-2]$yk™)dvX'  -  i 

Umx 


+  Atyy(X  +  2/x) 

At  X 

UCc 

f  C(|$|)(ry[x,  /- 2]$xk™  +  Yy[x,  l-2]Syk™)dVX'  -  1 


f  C(|$|)(ry[x,  /-I  \^xk'xy  +  yy[x,  l—\]^yklyy)dVX’  -  1 


+  A  tXyx 


■  (2-31) 


Finally,  (2-30)  can  be  discretized  spatially  using  a  simple  one-point  integration  and  point  match  testing, 
giving,  for  the  a -component  of  displacement. 


ux[Xi,  l ]  =  ~(yx  +  yy)ux[Xi,  l- 1]  -  YxYy“x[Xi,  1-2] 

At  Ni 

+  -E  CQbjMAxiJ-mfZ,  +  * T[x, ,  l-l]k%x)$xVj 

P  7  =  1 

At  Ni 

+  —  E  CQSijWAxj, l-\]kfxx  +  fT[xj,  l-\]kfyx)$xVj 

p  j= 1 
At  N> 

+  ~Yy  E  CQ&jMAxi,  l~2]k™x  +  \lrT[xit  l-2]k™  )$xVj 

P  7  =  1 

+  —  Yy  E  CQ&jWAxj,  l—2\kjVxx  +  ifT[Xj,  l-2]k™x)$xVj 

P  7  =  1 

At  N' 

+  —  E  C(|^|)(^[Xi,  l~l]k™y  +  fr[Xi,  l-l]kf^)$yVj 

P  7  =  1 

+  —  E  c(\^mx[xj,  i-mfiy  + 1 mxj,  i-\]kf;y^yVj 

p  j= i 

At  /v' 

+  —Yx  E  C(|^|)(^[Xi,  /— 2]/:- J  +  ^T[Xi,  l-2]k™  )$yVj 

P  7=1 

At  /v' 

+  —  Kx  E  CQ&jWAxj,  l~2]kfxy  +  iMx,-,  /— 2] k'J^y^yVj,  (2-32) 

P  7  =  1 

where  V)  is  the  volume  of  node  j,  Nj  is  the  number  of  nodes  in  the  neighborhood  of  node  i,  and 

$ij=Xj-Xi. 
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3.  Results 

The  PML  was  tested  on  two  types  of  problems,  first  a  wave  propagation  problem,  to  demonstrate  the 
effectiveness  of  the  PML,  and  second  a  crack  propagation  problem. 

3A,  Wave  propagation.  The  PML  was  first  tested  on  a  wave  propagation  problem  with  PML  boundary 
layers  and  a  Gaussian  distribution  as  an  initial  condition.  Specifically,  the  x -directed  displacement  was 
set  to 

ux(x,  t  =  0)  =  e-200|*-/>.nidl2,  (3-1) 

where  pmid  is  the  midpoint  of  the  region,  which  in  this  example  was  defined  as  0  <  x ,  y  <  1  and  discretized 
with  Ax  =  Ay  =  0.01  m.  The  Young’s  modulus  for  the  region  was  set  to  1  Pa,  the  Poisson's  ratio  was 
and  the  density  was  1  kg/m3.  The  PML  region  was  defined  as  the  0.3  m  border  around  the  lmx  1  m 
region  and  used  a  Gaussian  ramp  with  a  width  of  0.2  m,  finally  reaching  a  maximum  of  50  s"1  for  the 
remaining  0. 1  m.  For  the  Gaussian  kernel,  a  horizon  size  of  <5  =  1.1  Ax  was  used,  and  for  the  Heaviside 
kernel,  a  horizon  of  8  =  3.1  Ax  was  used.  The  kernel  constant  c  in  (2-4)  and  (2-5)  is  set  to  1  throughout. 

The  simulation  was  run  with  both  the  Heaviside  and  Gaussian  kernels,  with  the  total  strain  energy 
shown  in  Figure  2.  The  Gaussian  kernel  (the  dotted  line)  shows  the  largest  drop  in  energy,  reaching  a 
minimum  of  5.6  x  1 0  7,  and  the  Heaviside  kernel  (the  dashed  line)  decreases  to  1.  x  10~4.  A  bounded 
simulation  is  shown  for  reference  (the  solid  line),  which  used  a  fixed  displacement  boundary  condition 
and  the  Gaussian  kernel.  Figure  3  shows  a  waterfall  plot  of  the  x -directed  displacement  along  the  y  = 
0.5  line  for  the  Gaussian  kernel  with  the  PML  function  <px  shown  in  gray  on  the  far  end  of  the  plot 
(corresponding  to  t  =  1  s).  Figure  4  shows  the  absolute  value  of  the  x-directed  displacement  at  the  edge 
of  the  PML  region,  in  simulations  terminated  by  a  PML  and  with  a  fixed  boundary  condition.  The  wave 


Figure  2.  Total  strain  energy  in  a  simulation  terminated  by  a  PML. 
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<b  =  50  s-1 

X 


Figure  3.  The  x -directed  displacement  at  y  =  0.5  m,  terminated  by  a  PML. 


Figure  4.  The  x -directed  displacement  at  x  =  0.3,  y  =  0.5  m.  The  solid  line  shows 
results  terminated  by  a  PML,  and  the  dashed  line  used  a  fixed  boundary  condition. 


is  absorbed  at  the  boundary  with  minimal  reflections:  as  can  be  seen,  the  plots  align  for  a  time,  and  where 
they  deviate  (indicating  a  reflection  from  the  hard  boundary),  the  PML  simulation  remains  in  decay. 

For  verification,  the  method  was  compared  with  an  exact  analytical  solution.  Consider  a  cylindrically 
symmetric  wave  propagating  in  an  infinite  elastic  medium  with  the  same  constitutive  parameters  as  the 
above  example,  and  with  an  initial  condition  given  by 


(3-2) 
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Figure  5.  The  x-directed  displacement  at  x  =  0.5,  v  =  0  m.  The  solid  line  shows  the 
exact  solution  and  the  dashed  line  results  terminated  by  a  PML. 

where  here  we  take  b  =  1  and  a  =  0.1.  The  exact  solution  is  given  by  [Eringen  and  Suhubi  1975] 

Ur  , -  r- r2t2  I  4  r2t2 

u(r,  t)  —  —= - \/ R2  +  a( 2a  -  R2),  a  —  1  H - - - ,  R2  —  .  a2  -| - — ,  (3-3) 

sfl  aR6  a2  V  a 2 

where  c  is  the  longitudinal  wave  speed.  This  problem  was  simulated  in  a  two-dimensional  region,  2  m  x 
2  m  and  Ax  =  Av  =  0.01  m,  terminated  by  a  PML  with  the  same  dimensions  and  magnitude  as  the 
above  problem.  The  Gaussian  kernel  was  used  with  a  horizon  size  of  8  =  0.75  Ax,  with  an  actual  cutoff 
of  0.028  m.  The  results  are  shown  in  Figure  5,  with  the  exact  solution  shown  as  the  solid  line  and  the 
peridynamic  solution  shown  as  the  dashed  line.  The  peridynamic  solution  shows  good  agreement  with 
the  exact  solution  and  minimal  reflections  from  the  PML  boundary. 

3B.  Crack  propagation.  Crack  propagation  in  a  half-space  can  be  useful  for  modeling  physical  phenom¬ 
ena  such  as  indentation  experiments.  As  an  example,  we  model  such  a  problem  as  a  body  force  applied 
to  a  finite  region  with  small  precracks  in  a  region  terminated  on  three  sides  with  PMLs.  One  addition  to 
the  algorithm  for  this  problem  was  a  drag  term,  used  to  reduce  noise.  For  crack  problems  with  a  sudden 
force  application,  noise  and  oscillations  can  cause  hot  spots  and  undesirable  cracking.  To  remedy  this,  a 
drag  term  can  be  added  to  smooth  oscillations,  by  adjusting  the  nodal  velocity  as 

Nh 

v*[*<,  /]  =  (!-  D)v[Xj,  v[xj,  /],  (3-4) 

7  =  1 

where  D  is  the  drag  coefficient  and  /V/,  is  the  remaining  number  of  bonds  in  the  family  of  node  n  [Becker 
and  Lucas  2011]. 

An  absorbing  boundary  ensures  that  no  reflections  from  the  boundaries  interfere  with  the  crack  prop¬ 
agation,  possibly  causing  it  to  deviate.  Figure  6  gives  a  schematic  of  the  problem:  the  extent  of  the 
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Figure  6.  Schematic  diagram  of  the  problem  for  crack  propagation  in  a  half-space. 

computation  region  is  designated  by  the  solid  line,  the  PML  ramp  begins  at  the  dashed  line,  and  the  PML 
plateaus  at  the  dotted  line.  The  computational  region  was  70  mm  wide  and  35.25  mm  high,  the  PML 
region  began  at  15  mm  from  each  edge  (except  the  top)  and  peaked  at  5  mm  to  a  value  of  5  x  106s_1. 
The  node  spacing  was  0.496  mm  and  the  time  step  size  was  1  ns.  For  material  values,  the  density  was 
2235  kg/m3,  the  Young’s  modulus  was  65  GPa,  the  Poisson’s  ratio  was  0.2,  and  the  fracture  criteria  used 
a  fracture  energy  of  204  J/m2.  The  failure  criteria  used  in  this  simulation  was  bond-based,  that  is,  a  bond 
failed  if  it  was  stretched  past  a  given  limit,  determined  by  the  fracture  energy  [Ha  and  Bobaru  2010], 
The  maximum  relative  bond  stretch  was  then  2.971  x  1 0  3.  The  load  was  applied  across  a  10  mm  region, 
centered  at  the  top  surface,  with  precracks  on  each  edge  with  a  length  of  two  nodes  or  0.993  mm.  The 
simulation  was  run  for  a  total  of  10  /xs,  and  the  cracks  were  measured  manually  from  the  edge  of  the 
precracks  to  the  extent  of  the  damaged  area. 

Figure  7  shows  a  result  which  had  an  applied  load  of  250  N,  yielding  a  2.1 1  mm  crack.  Figure  8  shows 
a  close-up  of  the  damaged  area  from  Figure  7.  As  can  be  seen,  the  crack  extends  three  nodes  down  and 


Damage  (10000/10000)  Y 

0  0.217  0.433  \z  X 


Figure  7.  Damage  map  resulting  from  a  250  N  applied  load  after  10  /xs. 
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Figure  8.  Close-up  of  damage  map  from  Figure  7. 

three  nodes  across.  Finally,  the  applied  load  was  varied  between  140  N  and  500  N,  with  the  distance 
between  the  crack  tips  (the  crack  separation)  versus  the  applied  load  shown  as  the  dots  in  Figure  9. 
A  curve,  shown  as  the  solid  line  in  Figure  9,  was  fit  using  the  form 

d  =  Aps  +  B,  (3-5) 


Figure  9.  Crack  separation  versus  applied  force  for  indentation  into  an  elastic  half-space. 
Dots  represent  data  points  from  the  peridynamic  simulation  and  the  solid  line  is  a  curve 
fit. 
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where  d  is  the  crack  tip  separation  distance,  p  is  the  applied  load,  and  A  =  8.49  x  10~3  mm/N1,  5  =  1.16 
and  B  =  8.15  mm  were  determined  using  least  squares.  The  norm  of  the  residual  for  the  curve  fit  was 
1.65  x  itr6. 


4,  Conclusions 

A  perfectly  matched  layer  (PML)  was  applied  to  peridynamics  in  two  dimensions,  allowing  for  the  sim¬ 
ulation  of  infinite  regions.  State -based  peridynamics  was  used  as  more  flexible  constitutive  relations  arc 
necessary  to  implement  a  perfectly  matched  layer,  essentially  an  artificial  anisotropic  absorbing  material. 
Standard  discretization  techniques  were  used:  one-point  integration  and  point-matching  for  the  spatial 
discretization  and  forward  Euler  for  the  temporal  discretization.  Results  show  that  the  PML  absorbs 
incoming  waves  and  results  in  minimal  reflections  at  the  boundary  between  the  absorbing  layer  and  the 
computational  region.  A  Gaussian  function  was  used  as  a  ramp  to  avoid  these  numerical  reflections. 
Finally,  a  crack  propagation  problem  was  simulated  in  a  half-space,  modeling  indentation  problems. 
A  three-dimensional  implementation  of  the  method  would  be  straightforward,  following  an  identical 
procedure  for  the  peridynamics  and  PML  formulation. 
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